…recalculation

Major improvements to RTK positioning implementation:

- Fix filter initialization to use RINEX header position directly
- Implement ECEF to geodetic coordinate conversion (WGS84)
- Recalculate geometric ranges in updateFilter using current filter state
- Fix ambiguity initialization to use current geometric range
- Add outlier rejection for state updates > 10km to prevent divergence
- Add convert_rtklib_to_kml.py utility for visualization

Results:
- 118 stable FLOAT solutions (no divergence)
- 0.3-0.4m horizontal accuracy compared to RTKLIB reference
- Outlier rejection successfully prevents catastrophic divergence

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <noreply@anthropic.com>

…ncy)

Major milestone: libgnss++ is now a fully self-contained modern C++17 GNSS
positioning library with no runtime dependency on RTKLIB.

RTK Performance:
- Kinematic (1.2km): 100% fix rate, 12mm RMS (exceeds RTKLIB 98.3%)
- Short static (36m): 99% fix rate, 90mm RMS
- Long static (3.3km): 52% fix rate with WL-NL AR

New self-contained implementations:
- Eigen-based Kalman filter (algorithms/kalman.hpp)
- LAMBDA integer least-squares (algorithms/lambda.hpp + lambda.cpp)
- Saastamoinen troposphere model (models/troposphere.hpp)
- Klobuchar ionosphere model (models/ionosphere.hpp)
- Coordinate conversions with geodist (core/coordinates.hpp)
- Consolidated constants (core/constants.hpp)
- Solution writer (io/solution_writer.hpp)
- Native SPP with WLS + iono + trop (spp.cpp)

RTK algorithm improvements:
- Separate rover/base satellite positions from respective pseudoranges
- Analytical Sagnac correction in geodist()
- Wide-lane/Narrow-lane AR for long baseline ionosphere mitigation
- Melbourne-Wubbena cross-validation of LAMBDA integers
- Fix-and-hold with direct state constraint (PSD-preserving)
- Position-based hold validation
- High-variance DD pair exclusion (median-relative threshold)

RINEX 3 support:
- Navigation file parsing (G/R/E/C system chars, 4-digit year)
- Observation file parsing (epoch markers, per-satellite lines)
- Broadcast ionosphere parameter extraction (ION ALPHA/BETA, GPSA/GPSB)

Infrastructure:
- Regression test suite (tests/run_regression.sh)
- Analysis tools (tools/rtk_stats.py, plot_rtk.py, compare_rtklib.py)
- CLAUDE.md development guide
- Updated README.md with actual performance data

PreciseProducts::interpolateOrbitClock was returning false whenever the
query time matched a clock-only entry — the case for nearly every PPP
epoch when a 30-s IGS CLK file is loaded alongside a 15-min SP3 file.
The lower_bound short-circuit at line 870 ("if upper->time == time,
before=after") collapsed onto a CLK row whose position_valid was false,
so the function fell through to the single-sample branch and returned
false on the empty position. PPP's applyPreciseCorrections then silently
fell back to broadcast ephemeris for ~99% of receiver epochs.

Search the bracketing position-valid pair and the bracketing clock-valid
pair independently. Position is interpolated against the SP3 grid; clock
is interpolated against the CLK grid. New PPPTest regression covers a
query at an exact CLK-only timestamp.

Also iterate light-travel-time in PPPProcessor::applyPreciseCorrections
on the precise path so sat_position lies on the orbit at emission time
(broadcast path already iterated). The downstream geodist() applies the
Sagnac correction analytically.

Bench impact (TSKB DOY 105 / IGS final SP3+CLK, no ANTEX):
  baseline       3.23 m residual vs RINEX header
  this fix       3.82 m residual vs RINEX header
The +0.6 m drift exposes a separate downstream model bias (likely
clock-reference vs primary-signal mismatch); the broken precise path
was previously masking it via broadcast-eph fallback.

calculatePhaseWindup() existed in ppp_utils but was never called from the
PPP filter, leaving carrier-phase observables uncorrected for the
geometric rotation between satellite and receiver dipoles. The accumulating
~cm-class bias becomes substantial for sub-meter PPP and was contributing
to the residual 3-4 m offset measured at TSKB DOY 105.

Changes:
- calculatePhaseWindup now takes the sun position in ECEF, builds the
  satellite body frame from the spacecraft–Sun plane (yaw-steering
  attitude as used by IGS/RTKLIB), and uses an N–E–U receiver basis.
- The 2π ambiguity is resolved by rounding (previous − dphi) to the
  nearest integer rather than the prior `while (... > 0.5)` loop, which
  silently lost cycles when the inter-epoch step exceeded half a cycle.
- applyPreciseCorrections() computes the sun position once per epoch,
  then per satellite calls calculatePhaseWindup with the cached prior
  accumulator and subtracts the IF-combination correction
  windup_cycles × c/(f1+f2) from carrier_phase_if (or λ1 in
  per-frequency mode).

Bench: TSKB DOY 105 static PPP, IGS final SP3+CLK, RINEX header truth.
3D residual at tail: 3.82 m (PR #79 baseline) → 2.79 m (-1.03 m).

Adds two regression tests covering the 2π-ambiguity handling and the
zero-sun fallback.

The dual-frequency IF combination already cancels the first-order
ionospheric delay analytically (c1·I1 + c2·I2 = 0 by construction of
the IF coefficients). Subtracting an IONEX-derived correction on top
of an IF observable is mathematically a no-op but exposes the filter
to numerical-precision and TEC-mapping noise from the IONEX product.

Gate the IONEX subtraction in applyPreciseCorrections so it only fires
when (a) IF combination is disabled, (b) only a primary frequency is
available, or (c) per-satellite STEC estimation mode is on (which uses
IONEX as a tight state constraint, not an observation correction).

This is defensive — observed delta is small after PR for phase wind-up —
but cleans up the architecture so future iono-related model work
doesn't have to reason about a code path that should never fire.

Extend the existing receiver-only ANTEX loader so satellite blocks
(BLOCK IIA/IIR/IIF/IIIA/etc.) also load into PPPProcessor. Each entry
captures VALID FROM/UNTIL plus per-frequency PCO in the spacecraft body
frame (north→x, east→y, up→z; up = boresight toward Earth).

The matching application path is gated behind a new
ppp_config_.apply_satellite_antenna_pco flag, default OFF. Reason:
under the IGS convention switch in 2017, IGS final SP3+CLK products
are already published at the iono-free combination antenna phase
centre, so adding the body-frame PCO double-applies the offset.
Verified empirically on TSKB DOY 105 (2026-04-15): with the flag on
and the textbook (sat_pos += pco_ecef) sign, residual rises from
2.787 m → 3.207 m. A negated-sign run gives 2.633 m, but reproducing
the RTKLIB peph2pos sign — sign that DOES correct CoM-referenced SP3
to APC — empirically degrades the IGS final case, confirming the
convention finding.

Future SP3 sources known to ship CoM-referenced positions (some legacy
AC products, some GLONASS-only analyses) can opt in to the shift.

The loader itself is unconditional once --antex is supplied: it
populates satellite_antex_offsets_ for diagnostic/tooling use even when
the PCO is not applied to the geometry.

Investigation into the residual ~1.7 m gap to RTKLIB rnx2rtkp ppp-static
on TSKB DOY 105 (1.12 m vs our 2.79 m) localised the dominant
contributor: PPPProcessor::constrainStaticAnchorPosition() blends the
position state with a hard 50 % pull toward the SPP-derived anchor each
epoch and additionally zeros position×{clock, trop, ambiguity}
cross-covariances. This caps absolute static accuracy at the SPP
anchor's 2-3 m floor regardless of how clean the precise products are.

Disabling the blend on the precise-products path empirically diverges
(receiver-clock and zenith-trop states run away to 30 km within hours,
because the underlying observation model is leaking ~10–40 m systematic
into first-epoch pseudorange residuals).

This change just exposes the apply_static_anchor_blend knob (default
true, matching the legacy behaviour) and rewrites the comment so future
work has a clear place to start: identify the observation-model bias
that drives the divergence, then flip this flag false to lift the
anchor floor.

Bench unchanged at default (TSKB 3D = 2.787 m).

Two related observation-model fixes that close ~0.85 m of the gap to
RTKLIB rnx2rtkp on TSKB IGS final products (2.787 m → 1.934 m).

1. LTT (light-travel-time) extrapolation at SP3 boundaries
   PR #79 added a re-query of the precise-products interpolator at
   emission_time = reception − τ.  At the first epoch (reception time
   = first SP3 sample) the implied emission epoch falls 76 ms before
   the SP3 file starts.  The interpolator's single-sample fallback
   then returns sat_position at reception (no LTT correction) AND
   sat_velocity = 0, silently overwriting the bracket-derived
   velocity from the first call.  The geometric range was therefore
   off by sat_velocity·τ projected onto the line of sight — up to
   65 m on G01 at TSKB, which the EKF then absorbs into the trop
   state, walking trop to the 30 m clamp and forcing the receiver
   clock to swallow the rest.

   Replace the second SP3 query with a first-order Taylor expansion
   sat_position(em) ≈ sat_position(rcv) − sat_velocity·τ, which is
   accurate to mm for sub-second travel times and is well-defined at
   the SP3 boundary because sat_velocity from the first call is
   computed against the next bracket.  Per-sat std of the GPS prefit
   residual at TSKB epoch 1 collapses from 28 m to 8 m.

2. Phase-row gating threshold
   The phase rows in formMeasurementEquations were gated behind
   `convergence_min_epochs` (= 20 = 10 min at 30 s sampling) when
   precise products were loaded, vs `phase_measurement_min_lock_count`
   (= 1) on the broadcast path.  The static-anchor blend pins position
   to the SPP floor over those 20 epochs, so by the time phase rows
   activate the position covariance has been reset to 9 m² with all
   cross-covariances zeroed and the cm-class phase signal can no
   longer drive convergence.  Use the same threshold as the broadcast
   path.

Bench (TSKB DOY 105 IGS final SP3+CLK, 24 h static):
  baseline (PR #80 stack)       2.787 m
  + LTT Taylor extrapolation    2.065 m
  + phase from epoch 2          1.934 m
  RTKLIB rnx2rtkp ppp-static    1.121 m

GRAZ moves only marginally (2.20 → 2.19 m) because the LOS projection
of sat_velocity·τ depends on satellite geometry; TSKB's epoch-1 sky
view exposes the bug clearly.  Remaining 0.81 m gap to RTKLIB is
likely from the linear (vs Lagrange) SP3 interpolation away from
sample boundaries — separate work.

271 tests pass.

Replace the linear SP3 interpolator with a degree-9 Lagrange polynomial
(RTKLIB-style Neville's algorithm, 10-sample window) and disable the
static-anchor blend by default.

Linear interpolation between 15-min SP3 samples has theoretical sagitta
of 57 km at the chord midpoint (R·Δθ²/8 for GPS at Δθ = 7.5°); the LOS
projection produced km-class first-epoch pseudorange residuals (epoch 2
prefit: 4700 m on G01 at TSKB) that the static-anchor blend was
absorbing by pinning the position state at the SPP floor. With Lagrange
the same epoch-2 prefit drops to <11 m, the EKF converges cleanly with
no anchor needed, and the absolute residual approaches RTKLIB.

Bench (DOY 105 IGS final SP3+CLK, 24 h static):
                                  TSKB         GRAZ
  baseline (PR #80 stack)        2.787 m      2.200 m
  + LTT Taylor extrapolation     2.065 m      —
  + phase from epoch 2 (44c5f17) 1.934 m      2.186 m
  + Lagrange SP3 (this commit)   1.286 m      1.703 m
  RTKLIB rnx2rtkp ppp-static     1.121 m      —

The anchor flag stays in the config; it's still useful with broadcast
or SSR products whose orbit accuracy can be too poor for the filter to
converge unaided, but the default flips to false because Lagrange makes
the precise-products path RTKLIB-class without it.

Implementation notes:
- Neville's algorithm is numerically stable for the 8100 s × degree-9
  parameter range. Time origin is shifted to the query so the recursion
  evaluates at x = 0.
- Velocity is computed by re-evaluating the same polynomial at query+h
  (finite difference; degree-9 polynomial is smooth enough that h = 1 s
  matches the analytic derivative within the precision of the bench).
- The orbit and clock grids stay separate (SP3 at 15 min, CLK at 30 s);
  each interpolated by its own Lagrange window. The 30-s clock grid is
  smooth enough that the polynomial fit collapses to near-linear there.

271 tests pass.

Add VMF ATL multisite PPP bench fixtures

PPP: phase wind-up hook + IONEX gate + sat ANTEX loader